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ABSTRACT 

A nonlinear force-free solution is constructed for the coronal magnetic field in 
NOAA solar active region AR 10953 based on a photospheric vector magnetogram 
derived from Hinode satellite observations on 30 April 2007, taking into account 
uncertainties in the boundary data and using improved methods for merging 
multiple-instrument data. The solution demonstrates the 'self-consistency' pro- 
cedure of Wheatland & Regnier (2009), for the first time including uncertainties. 
The self-consistency procedure addresses the problem that photospheric vector 
magnetogram data are inconsistent with the force-free model, and in particular 
that the boundary conditions on vertical electric current density are over-specified 
and permit the construction of two different nonlinear force-free solutions. The 
procedure modifies the boundary conditions on current density during a sequence 
of cycles until the two nonlinear force-free solutions agree. It hence constructs 
an accurate single solution to the force-free model, with boundary values close, 
but not matched exactly, to the vector magnetogram data. The inclusion of un- 
certainties preserves the boundary conditions more closely at points with smaller 
uncertainties. The self-consistent solution obtained for active region AR 10953 is 
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significantly non-potential, with magnetic energy E/Eq ^ 1.08, where Eq is the 
energy of the reference potential (current-free) magnetic field. The self-consistent 
solution is shown to be robust against changes in the details of the construction 
of the two force-free models at each cycle. This suggests that reliable nonlinear 
force-free modeling of active regions is possible if uncertainties in vector magne- 
togram boundary data are included. 

Subject headings: Sun: corona — Sun: magnetic fields 



1. Introduction 



Intense magnetic fields around sunspot regions drive the most energetic examples of solar 
activity - large solar flares and coronal mass ejections - which can produce hazardous space 
weather conditions close to Earth. Severe space weather storms present a variety of hazards 

including increased radiation levels for space travellers and crew on polar air flights, the risks 

of dis abling of GPS systems and damage to power grids ( ICommittee On The Societal fc Economic Impacts 
20081 ). and potential large economic losses due to damage to communications satellites (Oden- 
wald et al. 2006). These effects motivate a need to understand and model the source regions 
of space weather threats at the Sun. 

Spectro-polarimetric measurements of magnetically sensi tive photospheric lines may be 
used to infer the vector magnetic field at the surface of the Sun (ILandi Degl'Innocenti fc Landolfi 
20041 ). and the resulting maps of the field are called vector magnetograms. A new genera- 
tion of instrumentation is set to provide vector magnetogram data over the next solar cycle 
with unprecedented quality, resolution, and temporal cadence, including the Solar Optical 
Telescope Spectro-Polarimeter (SOT/SP) on the Hinode spacecraft (Tsuneta et al. 2008) 
and the Helioseismic and Magnetic Imager ori the Solar Dynamics Observatory (SDO/HMI) 



( IScherrer. Hoeksema. &: the HMI Team! l2006l : iBorrero et al.l 120071 ) . In principle these mea- 



surements provide boundary values for coronal magnetic field modeling, or 'magnetic field 
reconstruction'. There is a pressing need for this modeling capability, which has many 
potential applications in solar research. However, studies have revealed basic difficulties 
preventing the construction of reliable coronal field models from vector magnetograms us- 
ing present modeling techniques. Specifically, a popular and well-developed approach to 
magnetic field reconstruction, nonlinear force-free modeling, has been shown to produce in- 



consistent results when apphed to Hinode/SOT vector magnetogram data (ISchrijver et al. 



i 



20081 : iDe Rosa et al.ll2009t IWheatland et al.ll2010[ ). Different modeling codes produce differ- 
ent results based on the same data, and the results from individual codes are also internally 
inconsistent. As described in more detail in section |2l vector magnetograms allow for two 
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sets of boundary conditions for the force-free model, and the two solutions do not in general 
agree. It is also difficult to construct the solutions: the boundary conditions on electric 
current present large currents which prevent the force-free methods from converging strictly. 
As a result the 'solutions' may be inaccurate with respect to the model. The modeling has 
been shown to be unable to provide a reliable estimate of the free magnetic energy of a solar 
active region. Further details of the problems are provided in section |2l 

One likely cause of the failure of nonlinear force-free modeling is a basic inconsistency 
between the boundary data and the model: the model includes only magnetic forces, an 
assumption that may be warranted in the corona, but which is thoug ht not to apply at 
the photospheric level where the data originate (iMetcalf et al.l Il995l ). Recently a new 
approach to coronal naagnet ic field reconstruction was presented to address this problem 
( jWheatland fc Regnierl |2009| ) . The 'self-consistency' procedure (described in more detail in 
section |2]) identifies a nonlinear force-free solution with boundary values close to those implied 
by the vector magnetogram data. The method constructs a single solution to the nonlinear 
force-free model which strictly converges. The solution does not ma tch exactly the field at 
the ph otospheric level, but may describe the coronal magnetic field. In lWheatland fc Regnier 
( I2OO9I ) the method was demonstrated in application to Hinode/SOT data for NOAA active 
region 10953. A solution to the model was constructed with energy E / Eq ^ 1.02, where Eq 
is the energy of the reference potential (current-free) magnetic field matching the boundary 
conditions on the vertical component of the field. The self-consistency procedure is capable 
of taking into account relative uncertainties in the solar boundary data, so that the bound- 
ary values for the solution match the vector magnetogram data more closely at points with 
small uncertainties. However, uncertainties were omitted from the IWheatland fc Regnier 
(120091 ) calculations, and this omission accounts in part for the low value of the magnetic free 
energy of the constructed solution (2% of the energy of the potential field), as explained in 
section 12.41 For this reason the results in IWheatland fc Regnien (120091 ) were argued to repre- 
sent a 'proof in concept' of the self-consistency procedure, rather than a full demonstration 
of its modeling capability. In this paper we return to the problem of modeling AR 10953, 
including uncertainties in the calculation. 

The detailed presentation of this paper is as follows. Section |5] provides a more detailed 
account of the background to the modeling presented here, incuding a summary of the 
nonlinear force- free model (section r2.1l) . a description of the failure of the model in application 
to solar data (section [2^2]) . and an account of the self-consistency proce dure (section [2l3|) 
and its initial application to AR 10953 by IWheatland &: Regnierl (120091 ) (section 12. 4p . A 
reader familiar with the background may skim or omit section [21 but it provides necessary 
background for a reader new to the topic. Section |3] describes the modeling of AR 10953 
taking into account uncertainties, with section 13.11 presenting the data, and section 13.21 the 
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results of the modeling. Section H] discusses the resuhs. 



2. Background 



2.1. Nonlinear force- free model 



Magnetic fields in the corona are often described as being force free ( 1Sturrocklll994j ) . The 



force-free model is static and includes only the magnetic (Lorentz) force, an approximation 



(Metcalf et al. 


1995; 


Gary 


2001) 



B ■ Va = and V x B = aB, 



'11 



where a is the force-free parameter. This parameter is constant along field lines, but varies 
(in general) from field line to field line, and the implicit dependence of a on B makes the 
model nonlinear. The boundary conditions required by the model are the normal component 
of B in the boundary (denoted -Bn), and a specification of a over one polarity (sign) of 
(e.g. Grad & Rubin 1958; Sakurai 1981; Aly 1989; Amari et al. 2006). We denote the 
two polarities P (positive) and (negative), for boundary points with i?„ > and i?n < 
respectively (assuming a normal directed into the coronal volume). The boundary condition 
on a is equivalent to a specification of the normal component of the electric current density 
J over either P or A^, because J = aB / /io in the model. 



(2007 


); 


Regnier 


(2007 


); 


Wieeelmann ( 


2008)1 



Neukirch 


(2005 


); 


Alv & Amari 



lem in a half space (2; > 0), with the lower boundary (2; = 0) representing the photosphere, 
solar curvature being neglected. The different methods vary in terms of numerical methods 
and in the treatment of boundary condition s. Some methods - e.g. the optimization method 



( jWheatland et al.ll2000l ; IWiegelmannl |2004| ) - use the full vector magnetic field over both P 
and A^ in the lower boundary as a boundary condition. This is formally an over-prescription, 
but does not introduce a problem provided the boundary conditions a re consistent with the 



model, as shown by application to test cases (ISchriiver et al. 



the 'current-field iteration' or Grad- Rubin procedure ( iGrad fc RubinI Il958t ISakurail Il981 



2006). Methods based on 



Amari et al.lll997l . |2006| ) use the formally correct boundary conditions on a descr ibed above. 
A Gr ad-Rubin implementation is the basis for the calculations in this paper (jWheatland 
20071 ). These methods work by replacing the nonlinear force- free equations by a sequence of 
iterations, such that at the i""^ iteration the linear equations 

1] 



B 



and V X B' 



(2) 
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are solved for given B^*"-^', subject to V ■ B^ = and the boundary conditions on 
and BW. For the half space problem the boundary conditions are the specification of at 
z = and a^^^ at z = over either P or A^. The iteration is usually started by setting B'^l 
equal to the reference potential (current-free) field matching the boundary conditions on Bz- 
A fixed point of the iteration scheme specified by Equations ([2]) provides a solution to the 
nonlinear force- free Equations ([1]). 



2.2. Failure of nonlinear force- free modeling for solar data 



Although the nonlinear force- free methods work well in application to test cases (jSchrijver et al 



20061 ) and to synthetic solar-like data (IMetcalf et al.ll2008l ). the application to vector mag- 
netogram data has revealed serious prob lems. Tests of differe nt methods usin g Hinode/SQT 



data for NOAA active regions AR 10930 (ISchriiver et al.ll2008r i and AR 10953 (IDe Rosa et al. 



20091 ) produced very different solutions for each active region depending on the specific 
method used. In particular the estimated magnetic free energy of each active region var- 
ied substantially among the solutions, preventing reliable determination of this quantity. 
Perhaps worse, the solutions from particular methods were also found to be internally in- 
consistent. For example, for the Grad- Rubin methods the vector magnetogram data provide 
both the P and choices for the boundary conditions on the force-free parameter. Specif- 
ically, a may be estimated from the data at the photosphere {z = 0) by calculating the 
gradients in the observationally-inferred transverse fields Bx and By according to 



a 



1 



dBy dBx 



dx 



dy 



(3) 



and these values are available over both P and A^. The P and solutions obtained by 
applying the Wheatland (2007) current-field iteration method to the data for active regions 



AR 1 0930 and AR 10953 were found to be very different (jSchrijver et al.ll2008l : iDe Rosa et al. 
2009h . 



It is likely that a number of factors contribute to the failure of the modeling. Vector mag- 
netogram field values are approximate due to resolution limitations, modeling assumptions, 
and photon noise in the polarization measurements. Assumptions made in the inference 



of the field from the Stokes spectral measurements (jPel Toro Iniesta fc Ruiz Cobd Il996 



Landi DegFInnocenti fc Landolfi 



introduce inaccuracy (lAuer et al 



2004 ) , for example using a Milne-Eddington atmos phere 



19771 : ISkumanich fc Liteslll987l : ISommier et al.ll2007h . The 



values are also uncertain due to the influence of instrumental noise and resolution limitations 
on the polarization measurements. The component of the magnetic field transverse to the 
line of sight is particularly uncertain due to the intrinsically smaller linear polarization signal 
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relati ve to the integrated intensity, magneto-optical effects (iLandolfi fc Landi Degl'Innocenti 
19821 ). and t he need to re s olve an intrinsic 180-degree anibiguit y in the azimuthal direction 
of the field jMetcaljllQoi iMetcalf et allbood : lleka et all 120091 ) . 



A problem with the force-free modeling is that different methods of solving the Equa- 
tions ([1]) make different assumptions regarding the problem of the 'missing' boundary con- 
ditions at the sides and top boundaries of the comput ational domain. The s e different as- 



sumptions lead (not surprisingly) to different results (ISchrijver et al.l l2006t iMetcalf et al. 



20081 : ISchrijver et al.ll2008t iDe Rosa et al.l 120091 ). However, the modeling also faces a more 



fundamental difficulty: the magnetic field at the height of the measurements is unlikely to 
be force free because the dense photospheric g as is subject to pres sure, gravity, and dynam- 
ical forces in addition to the magnetic forces (IMetcalf et al.lll995[ ). The boundary data are 
fundamentally inconsistent with the model. 

The accuracy of the force-free assumption for the boundary data may be checked a 
priori by calculating integrals of the field in the boun dary representing the net magneti c 



flux, the net force and the torque on the field in z > (iMolodenskiil Il969l : lAlyl Il984j . Il989l ) . 
The requirement that these integrals vanish provides a set of necessary conditions for the 
existence of a force-free solution for given boundary data. In general the integrals are found 
to be substantially different from zero for photospheric vector magnetogram data. One treat- 
ment of this problem is a procedure called 'preprocessing ' whereby the boundary data are 



varied to minimize the integra ls, subject to constraints (IWiegelmann. Inhester. &: Sakurai 



20061 : IWiegelmann et al.l 120081 ) . Preprocessing has been found to improve nonlinear force- 



free reconstructions from synthetic solar data for some methods (IMetcalf et al.ll2008l ). and 



preprocessed data were used in th e mo dehng of actiy e regio ns AR 10930 and AR 10953 



presented by ISchrijver et al.l (120081 ) and iDe Rosa et al.l (120091 ) respectively. However, there 



are non-trivial problems with the preprocessing approach. First, the conditions being min- 



(Aly 


1989; 


Reenier 


2007) 



and it is straightforw ard to show that prep rocessed boundary data are still inconsistent with 
the force- free model (IDe Rosa et al.ll2009l ). As a result, the P and solutions constructed 
by Grad-Rubin methods still disagree, and 'force-free solutions' constructed us ing other 
methods, e.g. the optimization method (IWheatland et al.ll2000l : IWiegelmannll2004l ) must be 
internally inconsistent with the force-free model. The situation is then that the observed 
boundary data is not matched by the 'solution,' and the force-free model is not being accu- 
rately solved. A second problem is that preprocessing typically involves a smoothing of the 
data, which is undesirable because it discards information. Recently, another solution was 
proposed, in which a vector field is calculated which attempts to simultaneously minimize the 
Lorentz force, the divergence, and the departure from observed boundary conditions, using 
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an optimization- like procedure ( iWiegelmann &: Inhesterll2010[ ). However, this method also 
fails to adhere to the force-free equations or retain (to within measurement uncertainties) 
the observed boundary data. 

A related problem with solar nonlinear force-free modeling is that the boundary con- 
ditions on the vertical electric current density corresponding to vector-magnetogram de- 
rived values of a include very large electric currents. The values of the current density 
are considerably uncertain because of the uncertainties in the field determinations discussed 
above, and there is also dependence of the values on the observational spatial resolution 



which changes w hen field values are interpolated onto new grids and rebinned ( lLekalll999 



Leka et al.ll2009l ). Nevertheless there appear to be electric currents sufficiently large to be 
inconsistent w ith the force-free model, which prevent str ict con vergence of n online ar force- 
free methods JSchriiver et al.l boosi : lOe Rosa et aD I2OO9I ) . The Iwheatlandl (12003) method 
is often found to reach a state where the field configuration exhibits slow oscillation as a 
function of Grad-Rubin iteration number rather than achieving a fixed state. The energy of 
the solution also oscillates slowly as a function of iteration number. This behavior suggests 
the non -existence, or the non-uniqueness, of a solution to the boundary value problem being 
solved (lAmari et al.l l2006l : llnhester fc Wiegelmanru |2006| ) . The consequence of this problem 
for the modeling is that the solution is non-unique (the exact result depends on the arbitrary 
choice of a stopping iteration) and also the force-free model is being solved less accurately 
(the field is less force free than in the case of strict convergence because the Grad-Rubin 
procedure does not achieve a fixed point). 



2.3. Self-consistency procedure 



In this paper we adopt the approach to resolving the inconsistency problem introduced 
recently by Iwheatland fc Regnierl ( 2009 ). The 'self-consistency' procedure identifies one 



nonlinear force-free solution with boundary values close, but not exactly matched, to the 
vector magnetogram data, with guidance from the estimated uncertainties. The result is an 
accurate solution within the assumptions of the force-free model. The procedure may be 
summarized as follows. 

First, P and nonlinear force- free solutions are constructe d, from unprepro cessed 
vector magnetogram boundary data, using the Grad-Rubin code of IWheatlandl (|2007| ). with 
a chosen number of iterations. These solutions use boundary values ± o"o for the force- 
free parameter, where the ao values are derived from the vector magnetogram data via 
Equation ([3]), and the uncertainties an are ob tained by propagating the uncertainties assigned 
to the field values (jLeka &: Skumanich 1999 ). The P solution provides a mapping, along field 
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lines, of the observational boundary values ao ± o"o at points in the region P to points in A^. 
This defines new possible boundary values ai ± o"i of the force-free parameter at points in 
N. (These values will not match the observational values ± ctq due to the inconsistency 
problem.) Similarly, the solution maps the observational boundary values ± ctq at 
points N to points in P, and so defines new possible boundary values ai ± (Xi at points in 
P. Together the two solutions define a complete new set of boundary values ai ± cri of the 
force-free parameter (over P and N). Figured] illustrates the first step of the self-consistency 
procedure. 



Second, Bayesian probability theory ( iJaynesI |2003[ ) is applied to decide on a boundary 



value of the force- free parameter at each boundary point, based on the two sets of values 
ao ± ctq and ai ± (Ti. The idea is that the values represent two observations subject to 
Gaussian uncertainties, and Bayes's theorem is applied to determine the true value of the 
force-free parameter. The resulting mo st probable value and associa ted uncertainty are given 



by the uncertainty- weighted averages ( IWheatland fc Regnieiil2009l ): 



The resulting set of boundary values 02 ± o"2 over P and is still, in general, inconsistent 
with the force-free model, but the values are expected to be closer to consistency. 

Finally, steps one and two are iterated. Two force-free solutions (the P and solutions) 
are constructed from the new boundary values 02 ± o"2, again by applying current-field 
iteration, in each case. The two solutions define new mappings between the P and A^ 
regions, and hence a complete new set of boundary conditions ± a^. Equations are 
applied to the two sets of values ^2 ± (^2 and ^3 ± to decide on most probable values 
04 ± CT4 at each boundary point, and so on. Each iteration of the self-consistency procedure 
[i.e. the construction of the P and A^ solutions, and then the application of Equations (jl]) 
to give a new set of boundary conditions on current] is referred to as a 'self-consistency 
cycle.' The cycles are numbered starting at one, and the P and A^ solutions are enumerated 
with an index k, such that the first cycle involves construction of P and A^ solutions for the 
vector magnetogram boundary values ao ±o"o, with the solutions numbered k = 1 and k = 2 
respectively. 

Self-consistency cycling is expected to converge to a self-consistent solution, i.e. a set of 
boundary values of the force-free parameter for which the P and A^ solutions are the same. 
The result is a single force-free solution, with (in some sense) a minimal departure from the 
observations. 
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2.4. Wheatland & Regnier (2009) 



In IWheatland &: Regnierl (|2009[ ) the self-consistency procedure was demonstrated in 
application to Hinode/SOT data for AR 10953 observed at 22:30 UT ori 30 Ap ril 2007. The 
vector magnetogram was the same as that described in iDe Rosa et al.l ( 20091 ) but without 
the preprocessing (and associated smoothing). Uncertainties were assumed to be equal at 
each boundary point. In that case the first of Equations (jl]) reduces to 02 = | («o + «i), i-e. 
a simple average of the two values. The self-consistency procedure was applied for ten cycles 
and found to converge to a solution with an energy E/Eq ^ 1.02, where Eq is the energy of 
the reference potential field. 



Figure |2]illustrates the self-consistent solution obtained by I Wheatland fc Regnierl (120091 ). 
Panel (a) shows the P solution and panel (b) shows the solution, after 10 self-consistency 
cycles. Each image looks down on the central part of the computational volume. The red- 
white-blue image in the background shows the boundary values of B^, with red indicating 
negative values, white values close to zero, and blue positive values. Two sets of field lines 
are shown in each panel by the red and the blue curves. The blue curves are field lines 
originating at points in the P polarity and the red curves are field lines originating at points 
in the A^ polarity. The sets of field lines in panel (a) a nd in panel (b) agree c l osely, illus- 
trating the achievement of self-consistency. Figure 1 in IWheatland fc Regnieii (120091 ) (not 
reproduced here) also shows the P and A^ solutions at the first self-consistency cycle, i.e. 
calculated directly from the original boundary data. The initial P and A^ solutions disagree 
markedly, with the A^ solution having more highly distorted magnetic field lines than the P 
solution, due to large electric currents associated with the leading polarity negative spot in 
the region. 

The self-consistent solution shown in Figure [2] has a relatively small rnagneti c free energy 
(about 2% of the energy of Eq). As discussed in IWheatland fc Regnieii (120091 ) . this is due 
in part to the neglect of uncertainties in the boundary values of a. Many of the boundary 
points have zero values of a. The Hinode field of view is small, and to compensate for this 
the Hinode data was embedded in a larger field of view based on a line-of-sight magnetogram 
from the Michelson Doppler Inter ferometer (MDI) inst rument on the Solar and Heliospheric 
Observatory (SOHO) spacecraft (jScherrer et al.lll995l ). The SOHO/MDI data points have 
zero values for and hence for a, and these zero boundary values are given equal weight 
to non-zero Hinode- derived values of a in the self-consistency averaging of the first of Equa- 
tions (jlj). This leads to a reduction in the electric currents in the solution, and hence a 
reduction in the magnetic free energy. Also, some strong currents are expected to occur in 
strong field regions, which will have more accurately-determined field values in the vector 
magnetograms (smaller uncertainties). These larger currents would be preserved in the ap- 
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plication of Equations (jl]) with un certainties included. This is e xpected also to contribute to 
a reduction in the currents in the IWheatland fc Regnierl (120091 ) solution. For these reasons 
the solution was presented as a proof of concept of the self-consistency procedure, rather 
than realistic attempt to model AR 10953. 



3. Modeling of AR 10953 including uncertainties 

3.1. Data 

The calculations presented here use a set of data for NOAA AR 10953 derived from the 
same Hinode/SOT SpectroPo larimeter (Tsuneta et al. 2 008) observations at 22:30 UT on 30 
April 2007 previously used in IPe Rosa et al.l and Iwheatland fc Regnierl J2009I ). The 

original Hinode data are ag ain inverted with t he HAO/ASP inversion code developed by Paul 
Seagraves at NCAR/HAO (IB. W. Litesll2009l ) which was used previously. The code provides 
estimates of the uncertainties in each magnetic component derived from the variation of 
the least-squares fit results for the observed Stokes spectra. These error estimates may be 
considered to provide a lower limit for the uncertainty in any single measurement. Standard 
error propagation is applied to the transformations used to obtain the field components 
x,By,Bz) in the Hinode/SP field of vi ew to give corresponding component uncertainties 



CTB^jCrByjCrB^) (ILeka &: Skumaniclj|l999l ). 



As with the iDe Rosa et al.l (120091 ) approach, co-aligned MDl line-of-sight data are used 
initially to produce a potential field with a very large field of view, including AR 10954 
roughly 15° to the east of AR 10953. However, the data-merging procedure deviates hence- 
forth. The two datasets are re-sampled onto a 0.5" grid (as compared to the original 0.32" 
and 1.8" resolution of the Hinode and MDl data, respectively). The computed MDI-derived 
image-plane potential vector field is merged with the coaligned Hinode/SP image-plane vec- 
tor field (each component separately), by including a Hanning-function weighting at the 
boundary between the two datasets. This approach ensures a smooth transition between 
the data over a few pixels and avoids introducing spurious vertical currents due to sharp 
gradients between the two datasets. 

The resulting combined Hinode/MDI vector field data and the uncertainties for each 
com ponent are cropped to a field of view similar in size to that used by lDe Rosa et al.l (120091 ) 



and Iwheatland fc Regnieii (120091 ) . and are re-sampled to the same 0.8" resolution used in 
the previous investigations, resulting in a final 313 x 313 grid. The values used here are 
fiux balanced to within 1%. 



Values of ao are calculated for all points with \Bz\ > 0.01 x ma,x{Bz) using Equa- 
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tion (|3]), with centered differencing used to approximate the derivatives. Uncertainties ctq for 
the force-free parameter boundary values are calculated at each point from (B^, By, By) and 



(ctr^ , aB ,0-3^ using propagation of errors applied to the differencing scheme (ILeka fc Skumanich 



19991 ). A nominal large uncertainty value, equal to the maximum over all points, is as- 
signed to points in the MDI field region. Points where the calculated ao is zero because 
l-B^I < 0.01 X ma.x{Bz) are also assigned this nominal large uncertainty. The effect of this 
is that poorly-determined boundary values of a carry little weight in the decisions on new 
boundary values made by the self-consistency procedure using Equation (jl]). 

Figure [3] illustrates the resulting vector magnetogram for NOAA AR 10953. Panel (a) 
shows the values of B^, and panel (b) shows the boundary values of J^. The same color table 
(red- white-blue for negative-ze ro-positive values) use d in Figure [2] is a dopted. Panel (b) may 



be com pared with Figure 3 in iDe Rosa et al.l (120091 ) and Figure 3 in IWheatland &: Regnier 



(120091 ) , which show the current density values for the vector magnetogram used in the earlier 
studies. The current structures in the two sets of boundary data are very similar, including 
strong regions of current in the leading negative polarity spot, and less concentrated current 
structures in the more diffuse region of positive polarity. The new current map lacks the 
artifacts associated with the edge of the Hinode data region visible in the earlier data. 

Figure H] illustrates the quantitative values of the uncertainty estimates for Bz and Jz 
obtained for the vector magnetogram. This figure displays non-parametric estimates (e.g. 



Silvermanlll986l ) of the density distribution of boundary points over values of \Bz\ and 
[panel (a)], and over values of |J^| and aj^ [panel (b)]. The uncertainty values for B^ 
illustrated in the panel (a) of Figure H] are relatively small overall because the Hinode/SP- 
derived uncertainties are based solely on fitting values, and do not include estimates 
of detection levels or systematic errors which would increase the uncertainty. Hence the 
estimates presented here should be considered to be lower bounds. The horizontal band of 
low uncertainty values for B^ in panel (a) of the figure represents the uniform noise in the 
MDI-derived areas. Panel (b) of Figure H] shows a substantial horizontal band of large values 
of \Jz\ with relatively small uncertainties, representing significant currents in the boundary 
data, as well as a vertical band of small \Jz\ values with large uncertainties, representing 
points where the current is poorly determined. 



3.2. Results 



The current-field iteration method (IWheatlandll2007l ) is applied to the data, as outlined 
in section 12.31 for 10 self-consistency cycles with 30 Grad- Rubin iterations used at each cycle 
to construct P and solutions. Each solution is calculated on a 313 x 313 x 300 grid. The 
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horizontal grid size matches the boundary data, and the vertical grid size is chosen to include 
all field lines which do not cross the sides of the computational volume. 

The treatment of 'missing information' associated with the absence of boundary condi- 
tions at the sides and to p of the coraputati onal volume is important for coronal magnetic 
field reconstruction. The IWheatlandl (120071 ) method addresses this problem by neglecting 
currents on field lines which cross these boundaries. At each iteration for each solution, 
points r = (x, z) on the grid threaded by field lines which cross the sides or top boundaries 
have a(r) set to zero. This includes points (x, y, 0) on the grid in the lower boundary, so this 
procedure modifies the boundary values of the force-free parameter during the Grad-Rubin 
iterations. For boundary points subject to this modification the corresponding uncertainty 
in the force-free parameter is set equal to the nominal large uncertainty value discussed in 
section 13. 1[ This choice minimizes the effect of the zero values during the construction [ac- 
cording to Equation (jlj)] of a new set of a boundary values at the end of a self-consistency 
cycle. The choice of the vertical grid size mentioned above also reduces the effect of this 
modification of the boundary conditions. 

The self-consistency procedure converges in less than ten cycles. Figure shows the 
self-consistent P and solutions [pan els (a) and (b) respectively] o btained after 10 cycles, 
using the same presentation as for the IWheatland fc Regnierl (120091 ) solution shown in Fig- 
ure |2l Once again the red-white-blue images in the background show the boundary values 
of Bz, and the blue and red curves are field lines originating in P and N polarities respec- 
tively. The sets of field lines shown in panels (a) and (b) agree closely, demonstrating the 
self-consistency achieved. Some specific differences are visible for field lines extending higher 
into the com putational volume. Ta ble [1] lists values of the vector field comparison metrics 
introduced by lSchrijver et al.l (120061 ) in the context of testing nonlinear force-free methods on 
known solutions. The more stringent of these checks, the mean vector error and the normal- 
ized vector error, are both about 2% for the solu tion compared with the P solution, which 



are small values for these met rics in this context (jSchrijver et al.l 120061 : iMetcalf et al.l 12008 



Wiegelmann fc Inhesterl |2010| ). The energies of the P and A^ solutions shown in Figure |5] 
are E/Eq ^ 1.0822 and E/Eq ^ 1.0819, which differ by less than 0.03%. The two fields are 
sufficiently similar that we refer to the result as representing a single self-consistent solution. 

The solution s hown in Figure |5] is considera bly more distorted from the potential con- 
figuration than the IWheatland fc Regnierl ( l2009l ) solution shown in Figure 121 consistent with 
the larger magnetic free energy (as discussed in section \2A\ the earlier solution has energy 
E/Eq fa 1.02). The energy of the self-consisten t solution obta ir ied he re {E/Eq ^ 1.08) is 
in the middle of the range of values reported in iDe Rosa et al.l ( 120091 ) for different nonlin- 
ear force-free methods applied to preprocessed boundary data (the values were in the range 
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E/Eq ^ 1.03-1.25). The energy of the self-consistent solution obtained here is larger than 
the energy E/Eq ^ 1.03 of the P solution constructed at the first cycle from the vector mag- 
netogram boundary values, but is smaller than the energy E/Eq ^ 1.17 of the solution 
constructed initially. The averaging of a values between foot points of field lines has led to a 
solution with an intermediate energy. The potential field has an energy Eq ^ 8.40 x 10^^ erg, 
so the magnetic free energy for the region is about 7 x 10^^ erg. 

Figure E] illustrates the convergence of the self-consistency procedure, showing the ener- 
gies of the P and solutions obtained at the end of each self-consistency cycle, versus cycle 
number. The P-solution energies are indicated by plus signs, and the iV-solution energies 
by diamonds. The energy of the first solution is substantially larger than the energy of 
the first P solution, as mentioned above. The solutions constructed at intermediate cycles 
approach one another in energy, converging in about eight cycles. During the initial self- 
consistency cycles, the Grad- Rubin iterations exhibit oscillation, as discussed in section 12.21 
The boundary data contains large values of the electric current density which prevent strict 
convergence of the Grad-Rubin method. This introduces some arbitrariness into the model- 
ing, because the exact P and N solutions, and hence the energies shown in Figure depend 
on the choice of a stopping iteration. This problem is discussed in more detail below. The 
Grad-Rubin solutions constructed for later cycles converge quite strictly (a fixed point of the 
iteration scheme is achieved). 

The self-consistency procedure changes the boundary values of the vertical electric cur- 
rent density and hence the boundary values of the horizontal field B/i = {B^, By), but does 
not change the boundary values of B^. Figure [7] illustrates the changes in the vertical cur- 
rent density Jz over the P and polarities in the boundary. The self-consistency procedure 
reduces the largest values of | | over both polarities, but retains some large currents. Many 
of the structures in the distribution of Jz observed in the initial data have survived, and 
in particular the patterns of strong-current in the leading negative polarity spot are pre- 
served [see panels (c) and (d) of Figure [7] , although the current distribution has changed 
in detail. Figure [8] illustrates the corresponding boundary values of the magnitude of the 
horizontal field at the photosphere, for the initial vector magnetogram and for the final self- 
consistent solution. Panel (a) of Figure |H] shows the horizontal field for the magnetogram 

(5*) + {By) , and panel (b) shows the the horizontal field for the self- 



data, Bl 



nl/2 



consistent solution, bI = {B^)^ + (-B^)' • The maximum value of the two fields is very 
similar, and there is a good qualitative correspondence between features in the two images, 
although the details have changed. 

Figure M illustrates the quantitative changes in the horizontal vector field, showing a 
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non-parametric estimate of the density of points (cf. Figure H]) over values of ABh/fJBt, as 

\2 I f-Df Di\2ll/2 



a function of the measured B^, where ^Bh. = [(-B^ — BD"^ + {B^ — -Bp^] and cr^^ is the 
uncertainty in B\^. Only points in the Hinode data region are included in this comparison. 
This figure indicates that the vast majority of points undergo small changes, but there is a 
substantial low-density tail of points subject to larger changes. The average absolute change 
in the horizontal field is (ABh) ^ 170Mx/cm^, and the average ratio of the changes to the 
uncertainties is (ABh/aB^) ~ 9. Many points are subject to changes in the horizontal field 
substantially larger than the corresponding uncertainty. As discussed in section 13. H the 
uncertainties assigned to the vector magnetogram data are lower bounds. Also, the initial 
boundary data are physically inconsistent with the force-free model (the departure from the 
model is not due only to observational uncertainties). The absolute changes in the horizontal 
field are comparable to those produced by preprocessing (jPe Rosa et al.ll2009l ). 



A measure of the inconsistency between the boundary data and the force-free model 
is given by the estimates of the net for ces and torques on the field i n z > which may 
be calculated from the boundary data (IMolodenskiil Il969l : lAlyl 11984 Il989l ). as discussed 
in section 12.21 Table [2] presents the values of the forces and torques for the initial vector 
magnetogram data, and for the self-consistent P and solutions. The components of the 
initial force are of order a few percent, and the components of the initial torque are of 
order 10%, in terms of relevant characteristic values. For the self-consistent solutions the 
force components are reduced to < 0.1%, and the components of the torques are reduced to 
< 0.1% (in most cases much less), with the same scaling. The changes in the horizontal field 
and hence current illustrated in Figures U\M are necessary to achieve this consistency with 
the force-free model. 

As mentioned above, the Grad-Rubin iterations do not converge exactly during the early 
and intermediate self-consistency cycles, and the choice of the final iteration for each cycle 
introduces some arbitrariness into the procedure. To investigate this effect the calculation 
is repeated, with A'gr = 20 Grad-Rubin iterations per cycle for both the P and A^ solutions 
(rather than A'gr = 30), also with A^gr = 40 iterations per cycle. The Grad-Rubin iterations 
converge in < 20 cycles for later cycles, so A^gr = 20 is approximately the minimum suitable 
choice. In each case, the self-consistency procedure is found to converge in less than 10 
cycles. The self-consistent fields obtained with A'gr = 20 and A'gr = 40 are quahtative very 
similar to the A"gr = 30 result, i.e. the field lines are visually very similar. The energies are 
also remarkably similar, with all three solutions having E/En ps: 1.08, to the s tated digits. 



Table [3] quotes the vector field comparison metrics from ISchrijver et al.l (|2006[ ) , comparing 
the P solutions of the two new cases (A'gr = 20 and A'gr = 40) with the original case 
(Agr = 30). The values are substantially larger than for the comparison of the individual 
P and A^ solutions (e.g. Tabled]), but are smaller than the typical differences introduced by 
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varying treatments of the unknown side and top boundary conditions in the application o f 



nonhnear force- free methods to known solutions (ISchrijver et al.ll2006l : iMetcalf et al.ll2008l ). 
Hence in the context of nonlinear force-free modeling, the solutions may be considered very 
similar. 

Figure [10] illustrates the boundary distributions of Jz for the three self-consistent 
solutions. Panel (a) shows the signal to noise ratio values Oq/'^o over the vector magne- 
togram, which determine how well the initial boundary conditions on are preserved at 
different locations in the boundary by the self-consistency procedure, according to the first 
of Equations (jl]). Panels (b), (c), and (d) show the boundary conditions on Jz in the self- 
conistent solutions for the cases A^^qr = 20, Nqb. = 30, and A^^qr = 40 respectively, using the 
same presentation as Figures E] and [71 Panels (b)-(d) are visually very similar. The same 
structures are observed, and the range of values of Jz is very similar. Comparison with the 
signal to noise ratio map [panel (a)] illustrates qualitatively how the current is unchanged at 
locations where ao (or Jz) is well-determined. These results suggest that the self-consistency 
procedure, including uncertainties, is identifying a particular force-free solution independent 
of this element of arbitrariness. 



4. Discussion 



A solution to the nonlinear force-free model is constructed for NOAA solar active re- 
gion AR 10953, based on a vector magnetogram derived from Hinode/SOT observations 
at 22:30 UT on 30 April 20 07. The solution applies the 'self-consistency' procedure of 
Wheatland fc Regnierl ( l2009l ). for the first time taking into account uncertainties in the 
boundary data. The boundary data are subject to improved techniques fo r data merging 



by co r nparison with earlier studies using the same Hinode/SOT observations (jPe Rosa et al. 



20091 : [Wheatland &: Regnierl |2009| ). The self-consistency procedure address es the problem 



that vector magnetogram da ta are inconsistent with the force-free model (jPe Rosa et al. 



2009[ : [Wheatland et al.[ [20101 ). In particular, the solar boundary data provide two possible 
force-free solutions (the P and A^ solutions, corresponding to choosing boundary conditions 
for the vertical electric current density Jz from the positive or the negative polarities of the 
boundary field respectively), which are generally found to be significantly different. The 
P and the A^ soutions are distinctly different for the vector magnetogram data used here. 
The self-consistency procedure modifies the boundary conditions on during a sequence of 
cycle s in which P and A^ solutions are constructed using Grad-Rubin iteration (jWheatland 
2007]), and arrives at a solution to the force-free model with boundary values close to, but 
not exactly matching, the observational vector magnetogram data. When uncertainties in 
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the boundary data are taken into account, the boundary values of Jz are preserved more 
closely at points having smaller uncertainties. 

Self-consistency modeling has two advantages over conventional force-free modeling. 
First, force-free methods tend to fail to strictly converge when applied to solar data, and 
hence do not accurately solve the force-free model. The self-consistency procedure achieves 
strict convergence. Second, the method identifies a single solution, rather than two solutions. 
Active region AR 109 5 3 was previously modelled using the self-consistency procedure in 
Wheatland fc Regnierl ( 120091 ) . but in that case uncertainties were not included, and the 
results were regarded as a proof of concept of the method, rather than a full demonstration of 
its capabilities. This paper presents a more comprehensive demonstration and investigation. 

The self-consistent solution obtained for active region AR 10953 is substantially non- 
potential, with a coronal magnetic energy E/Eq 1.08, where Eq is the energy of the 
potential field with the same vertical component of the magnetic field in the boundary. 
T his energy is significant l y larg er than the energy E/Eq k, 1.02 of the solution obtained 



m 



Wheatland fc Regnien (120091 ) . which is attributed to the neglect of uncertainties in the 



earlier paper. The non-potentiality of the region is due to strong currents in the inner parts 
of the region, in particular associated with the southern part of the the negative polarity 
leading spot, as shown in Figure 3. These currents were reduced in the earlier solution by the 
self-consistency averaging procedure [Equations (jl])] due to the omission of uncertainties (or 
rather, the treatment of all boundary points as having equal uncertainty). With uncertainties 
included correctly, these currents are preserved because they are in regions where the field, 
and hence current, is well-determined, as shown in Figure [TOl 

The energy of the new self-consistent solution is intermediate between the energies of 
the P and solutions constructed from the vector magnetogram data {E/Eq ^ 1.03 for 
the P solution, and E/Ec, ^ 1.17 fo r the N solution), and is also in the middle of the range 



of energies found by iDe Rosa et al.l (|2009[ ) using force-free codes apphed to boundary data 
derived from the same Hinode/SOT observations but subject to th e preprocessing procedure 
( jWiegelnaann. Inhester. fc Sakurail l2006l : IWiegelmann et al.ll2008l ). which is not used here. 
Recently ICanou &: Amaril (j2010[ ) presented a more detailed account of force-free m o deling 
of AR 10953 with two Grad-Rubin codes apphed to the A^-polarity iDe Rosa et al.l (|2009[ ) 
boundary data, and reported energies E/Eq ^ 1.27 and E/Eq ^ 1.31. The energies obtained 
by these authors are different from the energy of the initial A^ solution obtained here for a 
number of reasons. The Hinode/MDI boundary data used here is prepared differently, as 
discussed in section 13. H and preprocessing has not been apphed. The lack of convergence 
of the initial P and A^ solutions obtained here means that a range of energies is possible, 
depending on the choice of stoppping iteration (in practise we find energies for the A^ solution 
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in the range E/Eq 1.16 to E/Eq 1.22). The pr e sent m ethod treats the side and top 
boundaries in a different way to the ICanou &: Amaril (120101 ) method, and in particular we 
neglect currents on field lines w hich cross the side and to p bo undaries, which tends to reduce 
the energy. Finally, the stated ICanou fc Amaril ( 120101 ) and iDe Rosa et al.l (|2009[ ) energies 
refer to a smaller sub-region of the computational domain. Because of these differences we 
do not attempt a more detailed comparison with the earlier results. 

A degree of arbitrariness is introduced into the modeling by the lack of strict conver- 
gence of the Grad-Rubin iterations used to construct the P and solutions at early and 
intermediate self-consistency cycles. The results of the self-consistency procedure may de- 
pend on the choice of the number of Grad-Rubin iterations. This effect is investigated by 
repeating the calculation with A'gr = 20 and A'gr = 40 iterations at each cycle (rather 
than A'gr = 30). The two new self-consistent solutions are remarkably similar to the first 
solution, and in particular all three solutions have energy E / Eq ^ 1.08. This suggests that 
self-consistency modeling provides a solution to the problem of reliable estimation of the 
coronal magnetic free energy of an active region, once uncertainties in the boundary data 
are incorporated. The inclusion of uncertainties preserves boundary conditions on current 
at points in the boundary where the currents are most accurately determined. The three 
self-consistent solutions obtained with different numbers of Grad-Rubin iterations have very 
similar boundary distributions of Jz- 

The self-consistency procedure (including uncertainties) is a promising candidate for 
routine reconstruction of coronal magnetic fields. The technique is shown to produce a 
highly non-potential and accurately force-free coronal field model for NOAA AR 10953 from 
Hinode/SOT vector magnetogram data. The estimate of the free energy of the magnetic field 
appears to be robust based on the results obtained with varying numbers of Grad-Rubin it- 
erations per self-consistency cycle. The availability of high-quality vector magnetogram data 
over the next solar cycle from Hinode/SOT and SDO/HMI will provide ample opportuni- 
ties for further testing, development, and application of the method, to investigate many 
questions in the physics of solar activity. 
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Fig. 1. — The first step of tfie self-consistency procedure. Two solutions are constructed 
from the P (positive polarity) and N (negative polarity) boundary values of the force-free 
parameter ao, shown in panel (a). The P and N solutions map values of ao ± cq along 
magnetic field lines from P ^ N and from N ^ P respectively, as illustrated in panels (b) 
and (c) respectively. These mappings define a new set of boundary conditions (ai iui) over 
P and A^, as shown in panel (d). The blue curves in the boundary in each panel indicate 
the neutral line dividing P and A^, and the green curves in panels (b) and (c) are field lines 
of the P and solutions, with arrows indicating field direction. 



Fig. 2. — T he self-consistent solut i on for the coronal magnetic field of NOAA AR 10953 
obtained by IWheatland fc Regnierl (120091 ) . Panel (a) shows the P solution and panel (b) 
shows the N solution, after 10 self-consistency cycles. The color image in the background 
in each case shows the boundary values of S^, with negative values red, values close to zero 
white, and positive values blue. The blue curves in each panel are field lines originating at 
points in the P polarity, and the red curves are field lines originating at points in the N 
polarity. The sets of field lines for the P solution in panel (a) and for the N solution in panel 
(b) agree closely. 
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Fig. 3. — The vector magnetogram boundary data for NOAA AR 10953 used in the magnetic 
field reconstructions in this paper. Panel (a): photospheric values of the vertical field B^. 
Panel (b): initial photospheric values of the vertical electric current density J^. 
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Fig. 4. — The uncertainties in the vector magnetogram boundary values of the field, repre- 
sented using non-parametric estimates of the density distributions of boundary points. Panel 
(a): density distribution over values of \Bz\ and 0"^^. Panel (b): density distribution over 
\JJ and aj^. 



Fig. 5. — The self-consistent solution for the coronal magnetic field of NOAA AR 10953 
obtained using the new vector magnetogram boundary data including uncertainties. The 
presentation is the same as for Figure [2l Panel (a) shows the P solution and panel (b) shows 
the solution, after 10 self-consistency cycles. The color images in the background are 
the boundary values of B^, and the blue/red curves are field lines originating in the P/N 
polarities. The sets of field lines for the P solution in panel (a) and for the A^ solution in 
panel (b) agree closely. 
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Fig. 6. — The energy of the P solutions (plus signs) and the N solutions (diamonds) con- 
structed at each self-consistency cycle, in units of the reference potential field energy. 
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Fig. 7. — Initial boundary conditions on and the boundary conditions on for the self- 
consistent solution, for each polarity. Panel (a) shows the initial Jz values over P, and panel 
(b) shows the self-consistent Jz values over P. Panels (c) and (d) show the same quantities, 
respectively, over N. 
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Fig. 8. — Panel (a): The horizontal field at the photosphere for the vector magnetogram. 
Panel (b): the horizontal field at the photosphere for the self-consistent solution. 
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Fig. 9. — Quantitative changes in the horizontal field at the photosphere produced by 
the self-consistency procedure. A non-parametric estimate of the density distribution of 
boundary points, showing ABh/crB^ ^ function of the initial Bh magnitude, where 
aBf^ is the uncertainty in the magnitude of the initial horizontal field, and ABh = 
[{B^ — -B* )^ + {B^ — -By)^] is the change in the horizontal field, with i denoting values in 
the initial vector magnetogram, and / the final values in the self-consistent solution. 
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Fig. 10. — Panel (a) shows a map of the signal to noise ratio for ao, which determines how 
well the initial boundary conditions on are preserved by the self-consistency procedure. 
Panels (b), (c), and (d) show the boundary conditions on Jz for the self-consistent solutions 
obtained using 20, 30, and 40 Grad-Rubin iterations per self-consistency cycle respectively. 
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Table 1: Vector field comparison metrics for the calculated self-consistent P and N solution 
fields, for the default calculation (with 30 Grad- Rubin iterations per self-consistency cycle). 



Metric'^ Value 

1 - CS 8 X 10-6 

1 - VC 3 X 10-^ 

MVE 2 X 10-2 

NVE 2 X 10-2 



"These metrics were introduced by ISchriiver et all (|2006l ). CS is the Cauchy-Schartz metric, CS — 
Bii • 'B2i/ {BiiB2i), where Bi is the P solution and B2 is the N solution, the sum is over all grid 
points enumerated by i, and M = 313 x 313 x 300 is the total number of grid points. VC is the vector corrrela- 
tion, VC = J2^ Bh-B2, / yE^Bf^ET^- MVE is the mean vector error, MVE = Af"! J2^ \Bu-Bu\/Bu, 
and NVE is the normalized vector error, NVE ~ J2i l^u — B2i| /J2i ■ 



Table 2: Estimates of the components of the net force F and net torque T on the coronal 
volume above NOAA AR 10953. 



Quantity 


Vector magnetogram 


Self-consistent P solution 


Self-consistent N solution 




-1.2 X 10-2 


4.8 X 10"^ 


4.3 X 10"^ 


Fy/Fo 


3.9 X 10-2 


7.1 X 10"^ 


6.5 X 10"^ 


FjFo 


-7.0 X 10-2 


-5.1 X 10-^ 


7.1 X 10"^ 


r./Eo^ 


-1.2 X 10-1 


1.1 X 10-3 


1.9 X 10"^ 


r,/ro 


1.2 X 10-1 


-2.5 X 10"^ 


1.9 X 10"^ 


r./Eo 


3.9 X 10-2 


-7.6 X 10-^ 


-1.0 X 10-^ 



"The forces are in units of Fq = j pB dxdy with pB = (i?^ + By + i?^)/2/io, where the integral is performed 
over the magnetogram area with the origin for the coordinates at the lower left corner. 

''The torques are in units of Fq = / xpB dxdy, with pB — {B'^ + By+B^)/2fio, where the integral is performed 
over the magnetogram area with the origin for the coordinates at the lower left corner. 
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Table 3: Vector field comparison metrics for the self-consistent P solutions calculated with 
different numbers of Grad-Rubin iterations per self-consistency cycle. The comparisons are 
for the solutions with Nq^ = 20 and iVcR = 30, and for the solutions with A^^gr = 20 and 
A^GR = 30. 

Metric^ A^'gr = 30 and A^'gr = 20 A^'gr = 30 and A^'gr = 40 
1 - CS 5 X 10-^ 4 X 10-^ 

1 - VC 3 X 10-3 3 X 10-3 

MVE 0.08 0.09 

NVE 0.08 0.09 



"The definitions of the metrics are the same as given in Table [1] with Bi representing the A^gr — 30 solution 
and B2 the solution with either A^gr — 20 or A^gr = 30. 



